f 

COO-3077-145  { 


Courant  Institute  of 
Mathematical  Sciences 

ERDA  Mathematics  and  Computing  Laboratory 


A  Survey  of  Numerical  Methods 
for  Compressible  Fluids 

Gary  A.  Sod 


ERDA  Research  and  Development  Report 

Mathematics  and  Computing 
June  1977 


New  York  University 


UNCLASSIFIED 


ERDA  Mathematics  and  Computing  Laboratory 

Courant  Institute  of  Mathematical  Sciences 
New  York  University 


Mathematics    and   Computing  COO-3077-145 

A   SURVEY   OF   NUMERICAL  METHODS   FOR   COMPRESSIBLE   FLUIDS 

Gary  A.    Sod 
June   1977 


Contract    No.    EY-76-C-02-3077 


UNCLASSIFIED 


Abstract 

The  finite  difference  methods  of  Godunov,  Hyman,  Lax-Wendroff 
(two-step),  MacCormack,  Rusanov,  the  upwind  scheme,  the  hybrid 
scheme  of  Harten  and  Zwas,  the  antidiffusion  method  of  Boris 
and  Book,  and  the  artificial  compression  method  of  Harten  are 
compared  with  the  random  choice  method  known  as  Glimm's  method. 
The  methods  are  used  to  integrate  the  one-dimensional  equations 
of  gas  dynamics  for  an  inviscid  fluid.   The  results  are  compared 
and  demonstrate  that  Glimm's  method  has  several  advantages. 


I.  Introduction 

In  the  past  few  years  many  finite  difference  schemes  have 
been  used  for  solving  the  one-dimensional  equations  of  gas 
dynamics  for  an  inviscid  fluid.   Recently  the  random  choice  method 
(Glimm's  method)  introduced  by  Glimm  [6],  has  been  developed  for 
hydrodynamics  by  Chorin  [3].   Due  to  the  nonstandardness  of 
Glimm's  method,  as  well  as  the  difficulty  in  programming,  its 
acceptance  as  an  effective  and  efficient  numerical  tool  may  be 
restricted. 


In  the  following  sections  a  brief  discussion  of  the  methods 
is  given,  their  solution  to  a  sample  one-dimensional  problem  is 
compared,  the  advantages  of  Glimm's  method  are  discussed,  and 
finally  the  equations  used  by  Glimm's  method  are  derived  and  a 
flow  chart  for  the  programming  of  it  is  given. 

Basic  Equations.   The  one-dimensional  equations  of  gas  dynamics 
may  be  written  in  the  (conservation)  form: 

h^p  +^^(pu)  =  0  ,  (1) 

V+^x^^+P)  =  0  ,  (2) 


^m 
P 


V+^x^?  ^^+P))  =  °  (^^ 


where  p  is  the  density,  u  is  the  velocity,  m  =  pu  is  momentum, 

p  is  pressure,  and  e  is  energy  per  unit  volume.   We  may  write 

1   2 
e  =  pe  +2   pu  ,  where  e  is  the  internal  energy  per  unit  mass. 

Assume  the  gas  is  polytropic,  in  which  case 

P 


£  = 


r^Aip  '"^) 


where  ^  is  a  constant  greater  than  one.   Furthermore,  from  (4a) 
we  have 

P  =  A(S)p^  (4b) 

where  S  denotes  entropy. 

Equations  (l)-(3)  may  be  written  in  vector  form 

where 


u  = 


in 


1% 


and   F(U)  = 


m 
2 

H  (e+p) 


.; 


In  order  to  deal  with  solutions  which  contain  shocks,  we 
write  the  equations  in  integral  form,  which  is  obtained  by- 
integrating  equations  (l)-(3)  (or  equation  (5))  over  any  region 
in  the  upper  half  of  the  (x,  t)  plane  and  applying  Green's  theorem 
to  obtain  the  following  contour  integrals 


T  pdx  +  Cp  mdt  =  0  , 

f^  mdx  +  (f  (— +  p  )dt  =  0  , 

f  edx  +^  (-  (e+p))dt  =  0  . 


(6) 
(7) 
(8) 


II.  Description  of  the  Methods 
The  methods  of  Godunov  [5],  Lax-Wendroff  (two-step)  [I6], 
MacCormack  [18],  Rusanov  [polj  and  the  upwind  difference 
scheme  [19  ]  have  been  widely  used  and  no  benefit  can  be 
obtained  by  describing  them  here.   Hence,  these  schemes 
will  merely  be  listed  in  Table  I.   The  remaining  methods  under 
consideration  will  be  briefly  discussed. 


Glimm' s  Method.   Consider  the  nonlinear  system  of  equations  (5). 
Divide  time  into  intervals  of  length  At  and  let  Ax  be  the  spatial 
increment.   The  solution  is  to  be  evaluated  at  time  nAx,  where  n 
is  a  nonnegative  integer  at  the  spatial  increments  iAx, 
i  =  0,±1,±2,  ...  and  at  time  (n  +  -p)At  at  (i+-p)Ax. 


n 


The  method  is  a  two-step  method.   Let  u^  approximate 
U(iAx,nAt)  and  ^l^+Y/l   approximate  U(  (  i+1/2)  Ax,  (n+l/2)At)  . 
To  find  the  solution  u'?^M2  and  thus  define  the  method,  consider 
the  system  (5)  along  with  the  piecewise  constant  initial 
data 


U(x,nAt)  = 


n 


-i+1  ,    X  >  (i+l/2)Ax  , 


X  <  (i+l/2)Ax 


(9) 


Ax 
This  defines  a  sequence  of  Riemann  problems.   If  At<  2( I u] +c) ' 

where  c  is  the  local  sound  speed,  the  waves  generated  by  the 
different  Riemann  problems  will  not  interact.   Hence  the  sol- 
ution vCx,t)  to  the  Riemann  problem  can  be  combined  into 
a  single  exact  solution.   Let  g^   be  an  equidistributed 
random  variable  which  is  given  by  the  Lebesgue  measure  on 
the  interval  E-pjp"^-   Define 


n+1/2 
-i+1/2 


=  v((i+C^)Ax, (n+|)At) 


(10) 


At  each  time  step,  the  solution  is  approximated  by  a  piece- 
wise  constant  function.   The  solution  is  then  advanced  in  time 
exactly  and  the  new  values  are  sampled.   The  method  depends 
on  the  possibility  of  solving  the  Riemann  problem  exactly  and 
inexpensively . 


Chorin  [3]  (see  also  Sod  [21])  modified  an  iterative  method  due 
to  Godunov  [^]    which  will  be  described  below. 

Consider  the  system  (5)  with  the  initial  data 


s^  =  (pr'^rPp  '     X  <  0  , 

S^  =  (Pj,,u^,p^)  ,         X    >   0    . 


(11) 


The  solution  at  later  times  looks  like  (see  [1^])  Fig.  1,  where 

S-,  and  S„  are  either  a  shock  or  a  centered  rarefaction  wave. 

The  region  S^  is  a  steady  state.  The  lines  i      and  I      are  slip 

lines  separating  the  states.   The  slip  line  dx/dt  =  u^  separates 
the  state  S^  into  two  parts  with  possibly  different  values  of  p^, 
but  equal  values  of  u^  and  p^ . 

Using  this  iterative  method  we  first  evaluate  p^  in  the 
state  S^.   Define  the  quantity 

M   =  — .  (  12) 

^    u^  -  u^ 

If  the  left  wave  is  a  shock,  using  the  jump  condition  U,[p]  =  [pu], 
we  obtain 

M_g  =  p/^^-U^)  =  P*^^*-  U_g)  (13) 

where  U^  is  the  velocity  of  the  left  shock  and  p^  is  the  density 
in  the  portion  of  S^  adjoining  the  left  shock.   Similarly,  define 
the  quantity 

M  .^i-^.  (1*) 

r   u  -  u„ 
r   * 


If  the  right  wave  is  a  shock,  using  the  Jump  conditions 
U  [p]  =  [pu],  we  obtain 


M^  =  -p  fu  -U  )  =  -p,(u^-Uj 


(15) 


where  U  is  the  velocity  of  the  right  shock  and  p^  is  the  density 
in  the  portion  of  S^  adjoining  the  right  shock. 

In  either  of  the  two  cases  ((12)  or  (13)  for  M.  and  (l^) 
and  (15)  for  M  )  we  obtain 


M^  =  /pTpI  Mv^/vJ  , 


'r'^r 


(16a) 


where 


M„  =  /p.p.  ^{Vjv,) 


£^i 


y^   X  +  ^^   ,    X  >  1  , 


<l,(x)  =  < 


/-I 


1-x 


^^—^   ^  /o      >         X  <  1  . 

2/7  l-x-Y-l/^^ 


Upon  elimination  of  u^  from  (12)  and  (1^)  we  obtain 


(16b) 


(17) 


P*  = 


(u,  -  u  +— I  +  _l) 
"■    i        r  M^  M^  ^ 

TTx 


(18) 


Equations  ( l6a),  ( l6b),  and  (l8)  represent  three  equations  in  three 
unknowns  for  which  it  can  be  seen  that  there  exists  a  real  solution. 
Upon  choosing  a  starting  value  p^  (or  M  and  M  ),  we  iterate  using 

A/  -L 

equations  ( l6a),  (.l6b),  and  (l8).   For  details  of  the  starting 
values  see  Chorin  [3]  and  Sod  [2l]. 


After  p^,  M  ,  and  M  have  been  determined  we  may  obtain  u^ 
by  eliminating  p^  from  equations  (12  )  and  (1^), 

^* mttt^^ •  ^^^^ 

I        r 

The  finite  difference  method  due  to  Godunov  [5]  in  Table 
I  is  for  the  Eulerian  form  of  the  equations  of  gas  dynamics. 
The  method  developed  by  Godunov  [4]  for  the  Lagrangian  form 
is  also  a  two-step  method  where  the  second  step  is  the  second 
half  step  in  Table  I.   However,  the  values  of  li._|_-,  ^2  ^"^ 
p'?^^'^^  ^^^  replaced  by  \x^    (19)  and  p^^  (l8)  from  the  Riemann 
problem  at  i+1/2. 

Artificial  Viscosity.   In  the  methods  of  Godunov,  MacCormack, 
and  Lax-Wendroff  (two-step)  an  artificial  viscosity  term  was 
added.   The  artificial  viscosity  term  used  was  introduced  by 
Lapidus  [13].  It  has  the  advantage  that  it  is  very  easy  to  add  to 

an  existing  scheme  and  it  retains  the  high  order  accuracy  of  the 

'^n+ 1  /    \ 

scheme.   Let  u.    be  the  approximation  at  time  (n+IjAt  obtained  by 

any  one  of  the  above  schemes.   This  value  is  replaced  by  the  new 

approximation 

'^n   ~n  ~n 
where  A'u.  =  u.  -u._-|  and  v  is  an  adjustable  constant. 

This  equation  (20 )  is  a  fractional  step  for  the  numerical 

solution  of  the  following  diffusion  equation 

U   =  vAt(Ax5)[|u  (U  ]   . 
-t     Ax      "■  '  x'-x-'x 


It  is  shown  (see  Ladidus  [13]) that  this  new  difference  scheme 
(obtained  by  adding  the  artificial  viscosity)  satisfies  the  same 
conservation  law  that  the  previous  equation  did.   The  values  of 
the  constant  v  used  varied  from  method  to  method.   This  is 
discussed  in  the  section  on  numerical  results.  This  artificial 
viscosity  was  not  added  in  the  smooth  regions. 

Harten's  Corrective  Method  of  Artificial  Compression.   In  this 
section  we  discuss  the  Artificial  Compression  Method  (ACM) 
developed  by  Harten  [8],   This  method  is  designed  to  be  used  in 
conjunction  with  an  already  existing  finite  difference  scheme. 
The  purpose  of  this  method  is  to  sharpen  the  regions  which  contain 
discontinuities  whether  shocks  or  contact  discontinuities. 

Only  the  basic  idea  of  the  ACM  will  be  discussed  for  the 
case  of  a  single  conservation  law.   Let  u(x,  t)  be  a  solution  of 
the  conservation  law 

^t  "^^^^^x  =  °  ^^^  ^ 

which  contains  a  discontinuity  (u^(t),  u„(t),  S(t)),  where  u.  and 
Ut-,  are  the  values  on  the  left  and  right  of  the  jump  and  S  is  the 
speed  of  the  discontinuity.   The  discontinuity  is  either  a  shock 
or  a  contact.   Assume,  without  loss  of  generality  that  at  any 
given  time  t  the  solution  u  does  not  take  on  any  values  between 
UT(t)  and  u_(t).   Consider  the  function  g(u, t)  with  properties 

g(u,t)  sgn  [Uj^(t)  -u^(t)]  >  0  for  u  e  (u^(  t ),  Uj^(  t ) )  ,    (22) 
g(u,t)  =  0  for  u  ^  (u^(t),Uj^(t))  .    (23) 


This  function  g  will  be  called  an  artificial  compression  flux. 

It  can  be  seen  that  u  is  also  a  solution  of  the  conservation 
law 


u^  +  (f(u)+  g(u,t))^  =  0  . 


(2^) 


By  (23)  we  see  that  when  u  is  smooth  the  equation  (24)  is  identical 
with  equation  (21 )  and  the  shock  speed  S(t)  remains  the  same. 
Finally  it  is  observed  (from  (22))  that  if  (u  , u  , S)  is  a  shock 
or  contact  for  equation  (21 )  then  it  is  a  shock  for  the  modified 
equation  ( 24) . 

The  artificial  compression  method  solves  the  modified  equa- 
tion (24)  rather  than  the  original  equation  (21).  For  a  complete 
discussion  of  the  implementation  of  the  method  see  Harten  [8]. 

Let  u^    represent  the  approximate  solution  vector  to  (5) 
obtained  by  using  any  one  of  the  above  finite  difference  methods. 
In  solving  the  modified  system  (analogous  to  (2^  ) )  we  use  operator 
splitting.   We  first  define  the  difference  representation  g^  of 
the  artificial  compression  flux  g. 


li  =  "i^i 


(25) 


where 


and 


a.  =  max^  0,min 
k 


m 


in  (IS.+i/gUs^.i/g-Bgn  (5^^,/g))-|  ) 
1^1.1/2  1  ^1^-1/2  I        ^J 


(26) 


~   k       ~(k)  ~(k) 
where  k  refers  to  the  k-th  component  of  the  u,  ^■^^2_/2  ^   '^i+l~  ^i 
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Let  S.  -,  /p  represent  the  vector  whose  k-th  component  is  sgn  i^^^2_/2'^' 
Then  the  difference  scheme  which  applies  the  ACM  to  the  given  solu- 

~n+l  . ^ 
tion  u.    IS 

—X  - 

n+1       ~n+l       At    /^  ^        >, 

iii        =H         -2A3E   (gi+i-gi_i) 

+   2^    (l£l+l-Iil^i+l/2-  ISi-Si-ll^i.1/2)  (27a) 

'\'n+l  At    /„n  „n  -.  /' t7k  ^ 

=  ^1         -    2A7   (^i+1/2    -   ^i-1/2)'  (27b) 

,„n        nn      in      nio  t-^         <- 

where  9.1+1/2   "  %  "  %+l  "  'S-i+l  "  ^^  I  §-1+1/2'  applied  component- 
wise.  See  Harten  [7]- 

The  method  of  artifical  compression  is  designed  for  first 
order  schemes  and  cannot  be  applied  directly  to  higher  order 
schemes.   The  idea  of  ACM  is  based  on  the  existence  of  a  viscous 
profile.   See  Harten  [8].   Higher  order  schemes  introduce 
other  flux  terms  so  that  one  obtains  different  (nonphysical) 
speeds  of  propagation. 

Self-Adjusting  Hybrid  Schemes.   The  idea  of  self-adjusting  hybrid 

schemes  was  introduced  by  Harten  and  Zwas  [11].   Consider  a 

nonoscillatory  first  order  scheme  L-,  and  a  k-th  order  (k  >_  2) 

scheme  L,  , 
k' 


Vi  =   -1    -^^<.l/2-    ^-l/2\  (28) 

Vi  =  -i  -I7  (^1+1/2  -  ^Ll/2).  (29) 

So  as  not  to  violate  the  conservation,  hybridize  L,  and  L,  through 
their  numerical  fluxes  .   Define  the  hybrid  operator  L  by 


11 


L^  =  ^  -It  (f-i-Hl/2  -  ^-1/2^'  ^30) 

where 

^1  +  1/2  "  ^1+1/2  ^1  +  1/2  "^  ^1  -  ^1+1/2-*  ^1+1/2  '       ^31) 

6.,-,  /p  is  a  scalar  quantity  (called  a  switch)  which  satisfies 
0  ^  ^i+-|/p  5.  1  •   -^^  discontinuities  the  automatic  switch  is 
such  that  6  ~  1.   Hence  at  the  discontinuities  the  hybrid  scheme 
is  essentially  the  nonoscillatory  first  order  scheme. 
Equation  (30)  can  be  written  in  the  form 

L^  =  Vi  -^  17  ^^.1/2  (^1.1/2  -  ^1.1/2) 

-  ^-1/2  (^1-1/2  -  ^1-1/2^^  (32) 

so  that  if  6  is  o(Ax^)  where  the  solution  is  smooth,  then  for 
p  >_  k-1  we  have 

Lu.  =  L  u.  +  o(Ax^'^^)  .  (33) 

There  are  many  choices  for  such  schemes.   The  scheme  chosen 
here  is  discussed  in  Harten  [9] •   Taking  k  =  2  we  choose 
MacCormack's  scheme  and  by  adding  the  artificial  viscosity 
term 

i   ^^+1/2  (^1.1  -  ^i^  -  ^-1/2  (^1  -  ^1-1))         (3M 

to  MacCormack's  scheme  we  obtain  the  first  order  scheme. 
The  hybridized  scheme  becomes  for  the  system  (5) 


12 


uf' 


—1 


n  _  At  .pn    _  pH) 
-1    Ax  ^-i+1   -i^ 

At 


(35) 


1  (Uj^l  _  u?)  - 


2Ax  ^-i     -i-1^ 


-  }   <«l.l/2  (Itl.l  -  ^l)  -  «1-1/.  <"-! 


"-?-l" 


(36) 


The  stability  condition  for  the  first  order  scheme  is 

M  1^  ^  '^t  ^  /3 
max( |u I +c)  ^1  — ' 

this  being  stricter  than  the  stability  condition  for  MacCormack's 
scheme.   So  this  is  the  stability  condition  for  the  hybrid 
scheme . 

It  remains  to  describe  how  the  switch  6  is  chosen.   There 
are  many  possible  choices,  the  one  selected  is  described  in 
Harten  1 9J  •   Let  ^^  +  1/2   "  "^i+l  ~  ^i"   ^^^^"^ 


i+1/2 


-   A 


i-1/2 


^+1/2!  ' 


'i-1/2 


,  for  |A.^^/2 


,  otherwise. 


+  A 


i-1/2 


>  e 


(37) 


In  this  case  p  =  1  and  e  >  0  is  chosen  as  a  measure  of  negligible 
variation  in  the  density  p.   We  define  the  switch  9  by 


e 


1+1/2  =  "iaxce.,e.^^). 


Since  in  areas  which  contain  a  discontinuity  the  hybrid 
scheme  is  about  first  order  we  may  apply  the  artificial 
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compression  method  discussed  above.   However,  the  ACM  must 
not  be  used  in  smooth  regions.   For  this  purpose  the  switch 
is  used  again,  i.e.  equation  (27b)  may  be  replaced  with 


n+1    '\^n+l 


At  .  n     pn 
2Ax  ^"i+1/2  -i+1/2 


,n     pH     . 
1-1/2  -i-1/2^ 


(38) 


Antidiffusion  Method  of  Boris  and  Book.   In  this  section  we 
shall  discuss  briefly  the  antidiffusion  method  developed  by 
Boris  and  Book  [l].   The  purpose  of  this  special  technique 
known  as  "flux  correction"  is  to  achieve  high  resolution 
without  oscillations. 

It  can  be  shown  that  a  first  order  difference  scheme 
can  be  represented  by  an  equation  of  the  form 


^  +  fCu),  =  At  [gCu,f^)uJ^  , 


At. 


(39) 


where  gCu,^)  is  the  coefficient  of  the  diffusion  term. 

The  basis  of  the  antidiffusion  method  is  to  use  a  stable 
modification  of  a  diffusive  difference  scheme.   Let  the 
original  scheme  be  represented  by  (39),  the  modification  is 
represented  by 


u^  +  fCu)^  =  At 


/  /   At \     /   At  x  s 
(g(u,v-)  -  r(u,T^))u, 


'Ax' 


'Ax' 


(^0) 


where  r  is  a  positive  function.   One  can  introduce  the  anti- 
diffusion  term  by  operator  splitting.   The  first  step  consists 
of  solving 
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'vn+l     n 
with  the  original  difference  scheme,  say  u.    =  Lu..  Then 

In  the  second  step  let  A  be  a  difference  operator  approximating 
the  diffusion  equation 

-t  '   ^4^^"'i^  "x]x  =  0-  (^2) 

The  second  step  Is  the  antidlf fusion  step,  which  is  unstable 
by  Itself  since  It  approximates  the  backward  heat  equation. 
We  define 

n+1    n'^n+l    -T  n 
u.    =  Au.    =  ALu. . 

It  can  be  seen  that  If 

then  the  combined  scheme  AL  is  stable.   However,  (^3)  places 

more  of  a  restriction  on  7—  than  the  stability  condition  for  L. 

Ax 

We  chose  for  L  the  two-step  Lax-Wendroff  scheme.   Following 
Boris  and  Book  [  2]  ,  the  procedure  is 

u;:};M  (U?  ^  U^.,)  -^(F?,,  -Fj),         (..a) 

if  ^  =  5i  +  n(uj+i  -  2uJ  +  uj_^),  (44c) 

n+1  ^  ^n+1    .  c      _   c     s  (i^i^^s 

^1     -1      ^-1+1/2   -1-1/2^'  ^      ^ 

where 
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^i+l/2 


n+1   7?n+l. 


=  n  (u—  -u--). 


i+l/2 


=  U  Tt,  -  U..   , 


1+1 


—1 


^1+1/2  =  sgn(A.^^/2) 


max 


0, 


mm  |sgn(A.^^/2)^-l/2' 1^1+1/2! '"Sn(A.^^/2)A.^3/2]|. 

The  parameter  n  is  the  diffuslon/antldlf fusion  coefficient.   The 
stability  condition  is 

/ I   r    X  At 

max(  u  +c )  7 —  <  1 . 
'  '     Ax  — 

Hyman's  Predictor-Corrector  Method.   In  [12]  Hyman  describes 
a  predictor-corrector  type  scheme.   The  spatial  derivatives 
are  approximated  by  a  second  order  difference  operator  while 
the  time  derivative  (or  time  integrator)  uses  the  Improved 
Euler  scheme.   The  Improved  Euler  scheme  combines  a  first 
order  explicit  predictor  with  a  second  order  trapazoidal 
rule  corrector. 

For  stability  and  to  insure  proper  entropy  production  an 
artificial  viscosity  term  is  added.  The  artificial  viscosity 
term  used  Is  similar  to  that  used  by  Rusanov  [20]. 

The  scheme  is  given  by 


n+1/2 
u . 

—1 


=  u^  -  At  (DF^  -   6(*^^,/2 


4>"    )) 
"1-1/2^^' 


(i<5a) 


=  u .  -  At  P . 


n+1 

u. 

— 1 


=  u 


n  _  At  (DFfl/2  ^  pH)^ 


(^5b) 
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where 

o^  =  (u  +  c)^, 

and  c  is  the  local  sound  speed. 

The  stability  of  the  scheme  depends  on  the  number  of 
applications  of  the  corrector  (45b)  and  on  6 .   We  took  as  the 
stability  condition 

max  (|u|  +c)  ^  <  1. 
"      Ax   — 

In  order  to  maintain  stability,  the  artificial  viscosity  must 
not  be  completely  removed  in  the  smooth  regions.   However,  it 
can  be  reduced  in  these  regions  by  using  a  type  of  switch. 
The  one  chosen  was  suggested  by  Hyman  [12].   Replace  f^^+1/2    ^"^ 
(45a)  by  P  ^"+1/2  ^^^^^ 


P  = 


1   A  f   J^       >   n  ,  Ax 


1  ,  otherwise. 


This  type  of  switch  greatly  reduces  the  smearing  of  the 
contact  discontinuity  as  well  as  the  shock  wave.   This  switch 
is  a  type  of  artificial  compression. 

III.  The  Shock  Tube  Problem 

Figure  2  represents  the  initial  conditions  in  a  shock  tube. 
A  diaphragm  at  x^  separates  two  regions  (regions  1  and  5 )  which 
have  difference  density  and  pressures.   The  two  regions  are  in  a 
constant  state.   The  initial  conditions  are  p^  >  p   p-^  >  p  ,  and 
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u-|  =  u  =  0,  i.e.  both  fluids  are  initially  at  rest.   At  time 
t  >  0  (see  Fig.  3)  the  diaphragm  is  broken.   Consider  the  case 
before  any  wave  has  reached  the  left  or  right  boundary.   Points 
X-,  and  Xp  represent  the  location  of  the  head  and  tail  of  the 
rarefaction  wave  (moving  to  the  left).   Although  the  solution  is 
continuous  in  this  region  (region  2)  some  of  the  derivatives  of 
the  fluid  quantities  may  not  be  continuous.   The  point  x^  is  the 
position  that  an  element  of  fluid  initially  at  x„  has  reached  by 
time  t.   x^  is  called  a  contact  discontinuity.   It  is  seen  that 

across  a  contact  discontinuity  the  pressure  and  the  normal 
component  of  velocity  are  continuous.   However,  the  density  and 

the  tangential  component  of  velocity  are  not  continuous  across  a 
contact  discontinuity.   The  point  Xk  is  the  location  of  the  shock 
wave  (moving  to  the  right).   Across  a  shock  all  of  the  quantities 
(p,  m,    e,    and  p)  will  in  general  be  discontinuous. 

In  the  study  of  the  above  numerical  methods  the  following 
test  problem  was  considered:   p.,  =  1.,  p,  =  1.,  u-,  =  0.,  p   =  0.125, 
Pp-  =  0.1,  and  Uj-  =  0.   The  ratio  of  specific  heats  -y  was  chosen  to 
be  1.4.   In  all  of  the  calculations  Ax  =  0.01.   For  the 
Rusanov  scheme  the  value  of  w  was  taken  to  be  1.0.   In  the 
scheme  of  Boris  and  Book  the  parameter  n  was  taken  to  be  0.125. 
For  Hyman's  scheme  the  value  of  6  was  taken  to  be  0.8.   The  constant 
in  the  artificial  viscosity  term  v  was  taken  to  be  1.0  In  all 
but  one  case.   Also  the  value  of  o-  (see  Table  I)  was  taken  to  be  0.9- 

In  Glimm's  original  construction  a  new  value  of  C  was 
chosen  for  each  grid  point  i  and  each  time  level  n.   The 
practical  effect  of  such  a  choice  with  finite  Ax  is 
dlsasterous;  since  our  initial  data  is  not  close  to  constant 
Cwhich  was  an  assumption  made  by  Glimm).   In  fact,  if  C  is 
chosen  for  each  1  and  n,  it  is  possible  that  a  state  will 
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propagate  to  the  left  and  to  the  right  and  thus  create  a  spurious 
state.   An  improvement  due  to  Chorin  [3]  is  to  choose  |   only  once 
per  time  step  (hence  the  subscript  n).   The  details  of  the  method 
of  selection  of  the  random  number  are  found  in  Chorin  [3]  and 
Sod  [21] . 

Figure  4  indicates  the  results  using  the  first  order  accurate 
Godunov  scheme.   The  corners  at  the  endpoints  of  the  rarefaction 
wave  are  rounded.   The  constant  state  between  the  contact  dis- 
continuity and  the  shock  has  not  been  fully  realized.   The  transi- 
tion of  the  contact  discontinuity  occupies  7-8  zones  while  the 
transition  of  the  shock  occupies  5-6  zones. 

Figure  5  indicates  the  results  using  the  Godunov  scheme  with 
artificial  compression.   It  should  be  noted  that  for  this  case  the 
constant  in  the  artificial  viscosity  term  was  taken  to  be  2.0  to 
insure  that  the  solution  before  application  of  artificial  compres- 
sion was  oscillation  free.   For  the  artificial  compression  cannot 
be  applied  in  the  presence  of  oscillations.   The  corners  at  the 
endpoints  of  the  rarefaction  wave  are  still  rounded,  since  the 
artificial  compression  method  is  not  applied  in  smooth  regions. 
There  is  a  slight  undershoot  at  the  right  corner  of  the  rarefaction. 
Also  there  are  oscillations  at  the  contact  discontinuity.   The 
transition  of  the  contact  discontinuity  occupies  3-4  zones  while 
the  transition  of  the  shock  occupies  only  1-2  zones. 

Figure  6  shows  the  results  of  the  two-step  Lax  Wendroff  scheme. 
There  are  very  slight  overshoots  at  the  contact  discontinuity  and 
more  noticeable  overshoots  at  the  shock.   The  rarefaction  wave  is 
quite  accurate.   The  corners  at  the  endpoints  of  the  rarefaction 
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are  only  slightly  rounded.   The  transition  of  the  contact  dis- 
continuity occupies  6-8  zones  while  the  shock  wave  occupies  4-6 
zones.   It  is  observed  that  the  plots  in  Figure  6  are  quite  similar 
to  those  in  Figure  7  obtained  by  MacCormack' s  method. 

Figure  7  represents  the  results  of  the  second  order  MacCormack 
scheme.   There  are  slight  overshoots  at  the  contact  discontinuity 
and  more  noticeable  overshoots  at  the  shock  wave.   The  rarefaction 
wave  is  quite  accurate.   The  corners  at  the  endpoints  of  the  rare- 
faction are  only  slightly  rounded.   The  transition  of  the  contact 
discontinuity  occupies  7-8  zones  while  the  transition  of  the  shock 
occupies  5-6  zones. 

Figure  8  represents  the  first  order  accurate  Rusanov  scheme. 
The  contact  discontinuity  is  barely  visible  in  the  density  profile. 
The  corners  at  the  endpoints  of  the  rarefaction  wave  are  extremely 
rounded.   The  constant  state  between  the  contact  discontinuity  and 
the  shock  wave  is  barely  existent.   The  transition  of  the  contact 
discontinuity  occupies  l4-l6  zones  and  the  transition  of  the  shock 
occupies  6-8  zones.   This  scheme  is  extremely  diffusive.   This 
scheme  will  even  diffuse  entropy  for  zero  flow  fields. 

Figure  9  represents  the  Rusanov  scheme  with  artificial  com- 
pression.  The  results  with  artificial  compression  are  greatly 
improved.   The  corners  at  the  endpoints  of  the  rarefaction  wave  are 
still  rounded  since  the  artificial  compression  method  is  not  applied 
in  this  area.   The  constant  state  between  the  contact  discontinuity 
and  the  shock  is  much  more  visible.   The  transition  of  the  contact 
discontinuity  occupies  2-3  zones  while  that  of  the  shock  wave 
occupies  only  1-2  zones. 
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Figure  10  represents  the  upwind  difference  scheme.   It 
Is  observed  that  between  the  left  constant  state  and  the 
left  endpoint  of  the  rarefaction  wave  Is  a  shock  (discontinuity) • 
This  Is  clearly  a  nonphyslcal  solution.   This  Is  a  result  of 
the  method  used  to  staballze  the  scheme,  by  using  centered 
differences  for  the  pressure  term  In  the  momentum  equation. 

Figure  11  shows  the  results  of  the  Gllmm  scheme.   The  shock 
wave  and  the  contact  discontinuity  have  been  computed  with  Infinite 
resolution,  i.e.  the  number  of  zones  over  which  the  variation 
occurs  is  zero.   Due  to  the  randomness  of  the  method  the  positions 
of  the  shock  and  the  contact  discontinuity  are  not  exact.   However, 
on  the  average  their  positions  are  exact.   The  corners  at  the 
endpolnts  of  the  rarefaction  wave  are  perfectly  sharp.   It  is 
observed  that  the  rarefaction  is  not  smooth,  yet  it  is  extremely 
close  to  the  exact  solution.   The  constant  states  are  perfectly 
realized. 

The  Gllmm  scheme  requires  between  2  and  3  times  as  much 
time  (see  below)  as  the  other  finite  difference  schemes  tested. 
However,  the  Glimm  scheme  requires  far  less  spatial  grid  points 
for  the  same  resolution.   This  is  displayed  in  Table  II,  where 
9  Interior  grid  points  are  used.   All  details  are  visible. 

The  Gllmm  scheme  on  the  average  is  conservative.   One 
other  check  on  the  accuracy  is  to  use  the  conservation  laws 
(mass,  momentum,  and  energy).   For  example,  the  total  mass  is 
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evaluated  by 

Qp  =  2^D(iAx)  AX. 
1 

In  Table  III  the  values  of  the  total  mass,  momentum,  and 
energy  are  displayed.   The  mass  and  the  energy  are  seen  to 
be  conserved  on  the  average,  i.e.  there  are  fluctuations 
but  they  are  contained  within  a  small  interval.   The  momentum 
is  seen  to  increase  linearly  on  the  average  (allowing  for 
fluctuations ) . 

Figure  12  shows  the  results  of  the  antidif fusion  method 
of  Boris  and  Book  applied  to  the  two-step  Lax-Wendroff  scheme. 
There  is  a  slight  overshoot  at  the  right  corner  of  the 
rarefaction  .   The  rarefaction  wave  is  very  accurately  computed. 
The  corners  at  the  endpoints  of  the  rarefaction  are  only  slightly 
rounded.   The  constant  state  between  the  contact  discontinuity 
and  the  shock  wave  is  only  partially  realized.   The  transition 
of  the  contact  discontinuity  occupies  5-7  zones  and  the 
transition  of  the  shock  occupies  1-2  zones.   The  resolution  is 
much  better  than  the  two-step  Lax-Wendroff  scheme  alone 
(see  Fig .  6 ) . 

Figure  13  repsents  the  hybrid  scheme  (35)  and  (36)  of 
Harten  and  Zwas .   The  solution  is  free  of  oscillations.   The 
corners  at  the  endpoints  of  the  rarefaction  wave  are  only 
slightly  rounded.   The  constant  state  between  the  contact 
discontinuity  and  the  shock  is  only  partly  realized.  The  trans- 
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ition  of  the  contact  discontinuity  occupies  8-9  zones  and  the 
transition  of  the  shock  occupies  5-6  zones. 

Figure  l4  represents  the  hybrid  scheme   of  Harten  and  Zwas 
with  the  use  of  artifical  compression.   Since  the  artlfical 
compression  is  not  applied  in  smooth  regions  the  rarefaction 
is  the  same  as  in  Fig.  13-   The  transition  of  the  contact 
discontinuity   occupies  3-4  zones  and  the  transition  of  the 

shock  wave  occupies  2-3  zones. 

Figure  15  represents  the  results  of  Hyman's  predictor- 
corrector  scheme,  where  the  corrector  has  been  applied  once. 
The  solution  is  oscillation  free.   The  corners  at  the  endpoints 
of  the  rarefaction  are  almost  perfectly  sharp.   The  constant 
states  between  the  rarefaction  and  the  contact  discontinuity 
and  between  the  contact  discontinuity  are  extremely  well 
defined.   The  transition  of  the  contact  discontinuity  occupies 
6-8  zones  while  the  transition  of  the  shock  occupies  3-^  zones. 

The  timing   results  for  all  of  the  methods  are  listed  In 
Table  IV.   The  times  are  for  100  spatial  grid  points.   The  only 
substantial  'difference  in  timing  is  between  Glimm's  scheme 
and  the  other  finite  difference  schemes.   For  Glimm's  scheme 
requires  between  2  and  3  times  as  much  time.   However,  Glimm's 
scheme  can  give  the  same  resolution  with  far  less  points 
Cas  seen  in  Table  II).   From  the  point  of  view  of  the  least 
number  of  grid  points  per  desired  resolution,  the  Gllmm  scheme 
can  be  seen  to  be  much  faster. 

IV.  Conclusions 

Of  all  the  finite  difference  schemes  tested,  without  the 
use  of  corrective  procedures,  Godunov's  and  Hyman's  methods 
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produced  the  best  results. 

It  is  obvious  from  the  figures  that  the  Gllmm  scheme  gives  the 

best  resolution  of  the  shocks  and  contact  discontinuities. 
Glimm's  scheme  is  at  best  first  order  accurate  (see 

Chorin  [3])  so  that  boundary  conditions  are  easily   handled. 

It  is  possible  that  the  rarefaction  wave  obtained  by 
Glimm's  method  can  be   smoothed  out  by  a  type  of  averaging. 
This  is  presently  being  considered. 

The  hybrid  method  of  Harten  and  Zwas  combines  first  and 
high  order  schemes  in  such  a  way  as  to  extract  the  best 
features  of  both.   The  high  order  scheme  produces  better 
approximations  to  the  smooth  parts  of  the  flow. 

The  corrective  procedures  of  Boris  and  Book  and  Harten 
improve  the  resolution  of  a  given  scheme.   The  artifical  com- 
pression method  being  restricted  to  first  order  schemes  except 
when  used  in  conjunction   with  the  hyblrd  type  schemes 
produces  far  better  results  than  the  antldlf fusion  method  of 
Boris  and  Book.   Both  methods  are  easily  added  to  existing 
programs  (as  a  subroutine) .   The  antidiff uslon  method  requires 
slightly  more  storage  than  the  artificial  compression  method 
since  the  former  must  retain  two  time  levels  of  information 
for  the  computation  of  Intermediate  results  (equation  (44c)). 

A  major  disadvantage  of  the  antldlf fusion  method  of 
Boris  and  Book,  the  hybrid  scheme  of  Harten  and  Zwas,  and 
the  artificial  compression  method  of  Harten  is  that  there 
are  a  number  of  parameters  to  be  chosen,  which  depend  on  the 
given  problem.  In  the  antldlf fusion  method  the  coefficient 
of  diffusion/'antidiffuslon  must  be  chosen.   The  value  of  this 
parameter  can  greatly  affect  the  results.   In  the  hybrid 


24 


scheme  a  tolerance  must  be  chosen  for  the  automatic  switch  which 
is  taken  to  be  a  measure  of  negligible  variation  in  entropy  or 
density  for  example.   This  tolerance  depends  on  the  given  problem. 
In  the  artificial  compression  method  a  test  must  be  included  to 
locate  the  rarefaction  (and  other  smooth  regions).   Many  of  the 
standard  tests  fail  to  work  well  enough  for  the  use  of  artificial 
compression. 

With  the  method  described  for  solving  the  Riemann  problem 
in  the  Glimm  scheme,  it  can  only  be  used  for  the  equations  of  gas 
dynamics  in  rectangular  coordinates.   It  is  possible  to  generalize 
Glimm' s  method  to  other  coordinate  systems  and  different  equa- 
tions.  See  Harten  and  Sod  [10]. 

The  applicability  of  Glimm 's  method  to  other  geometries  has 
only  just  started  to  be  explored.   One  successful  application  is 
to  the  equations  of  gas  dynamics  for  a  cylindrically  or  spheri- 
cally symmetric  flow.   See  Sod  [  22]  . 
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Appendix:   Implementation  of  Glimm's  Method 

In  this  appendix  we  discuss  the  equations  required  for  the 
computer  implementation  of  Glimm's  method. 

As  in  Fig.  1,  the  fluid  initially  at  x  <  0  is  separated 

dx 
from  the  fluid  initially  at  x  >  0  by  a  slip  line  -grp  =  u^ .   There 

are  a  total  of  10  cases  to  consider, 

I.     The  sample  point  i   Ax  lies  to  the  left  of  the  slip  line 


'n 


(^^Ax  <  u^At/2). 

(a)  If  the  left  wave  is  a  shock  wave  (p^  >  p  )  and  (l)  if 

dx 
^  Ax  lies  to  the  left  of  the  shockline  -r-p  =  U  ,  we  have  p  =  p., 

u  =  u.,  and  p  =  p.,  (2)  if  ^  Ax  lies  to  the  right  of  the  shockline 

O  lb  i  i 

dx 

-ry  =  U,,  we  have  p  =  p^,  u  =  u^,  p  =  p^,  where  p^  can  be  obtained 

from  (13) 

M 


I 


U^-u. 


(^6) 


(b)  If  the  left  wave  is  a  rarefaction  wave  (p^  Jl  ?»)• 

Define  the  soiind  speed  to  be  c  =  j^.   The  rarefaction  wave  is 

dx 
bounded  on  the  left  by  the  line  defined  t»y  -^  =  u  -  c   ,    where 

i7Pfl  dx 

c.  =  — -,    and  on  the  right  by  the  line  defined  by  -^-p  =  u^  -  c^, 

where  c^  =  / .   The  flow  is  adiabatic  in  smooth  regions,  so  m 

*   '  P* 

this  region  A(s)  in  (4b)  is  a  constant,  denoted  by  A,  and  we 
obtain  the  isentropic  law  p  =  Ap^.   p^  is  obtained  by  using  the 
isentropic  law 

P^P^"^  =  P*P*^  =  A  .  (^7) 
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Then  we  obtain  from  {^Q 


(1)  If  I  Ax  lies  to  the  left  of  the  rarefaction  wave,  then 


'n 
p  =  p^,  u  =  u^,  and  p  =  p^. 

(2)  If  ^  Ax  lies  inside  the  left  rarefaction  wave,  we  equate 

dx 
the  slope  of  the  characteristic  -^^  =  u  -  c  to  the  slope  of  the  line 

through  the  origin  and  (^^Ax,At/2),  obtaining 


2i   Ax 
u  -  c  = 


""  {^9) 


At 
With  the  constancy  of  the  Riemann  invariant 

2c(7-l)-^  +  u  =  2c^(7-l)'^+u^  ,  (50) 

the  isentropic  law,  and  the  definition  of  c,  we  can  obtain  p,  u, 
and  p.   Using  the  isentropic  law  we  obtain 

p  =  p,p7p^  =  Ap^  •  ^^') 

Using  equation  ('50)  we  obtain,  by  solving  for  c 

c  =  c^  +  ^^  (u^-  u)  .  (52) 

By  substitution  of  (52)  into  {^9)    and  solving  for  u  we  obtain 

o    2P  Ax      /       1  \ 


u 

7" 


By  substitution  of  ( 53)  into  ( 52)   c  is  obtained;  by  substitution 
of  (  52)  into  the  definition  of  c  and  solving  for  p  we  obtain 
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P  =  (f^)'^"'"'  •  (53) 

(3)  If  Ij^Ax  lies  to  the  right  of  the  left  rarefaction  wave 
we  obtain  p  =  p^,  u  =  u^,  and  p  =  p^. 

II.    The  sample  point  ^  Ax  lies  to  the  right  of  the  slip  line 
(I^Ax  >  u^At/2). 

(a)  If  the  right  wave  is  a  shock  wave  (p^  >  p  )  and  (l)  if 

dy 
i   Ax  lies  to  the  left  of  the  shockline  defined  by  -^   =  U  ,  we  have 
n  dt     r 

p  =  p#j  u  =  u^,  and  p  =  p^,  where  p^  is  obtained  from  (15) 

-M 


^  r_ 

'*    u^  -  U 


(5M 


'■* 


r 


(2)  If  ^  Ax  lies  to  the  right  of  the  shockline  defined  by 
■^  =  U^,  we  have  p  =  p^,  u  =  u^,  and  p  =  p^, 

(b)  If  the  right  wave  is  a  rarefaction  wave  (p^  "^  P  ) .   The 
rarefaction  wave  is  bounded  on  the  left  by  the  line  defined  by 


dx  i^P* 

-TT-  =  u^  +  c^,  where  c^  =  j and  p^  can  be  obtained  from  the 

isentropic  law 

P„p"^  =  P*p;^  =  A  .  (55  ) 


Then  we  obtain  from  (55  ) 


P*  =  (-^)    ;  (56  ) 


dx  /^Pr 

and  on  the  right  by  the  line  defined  by  -rrr   =  u  +  c  ,  c   =  / , 

tz>  J  "'dtrr'rVp 

(1)  If  i   Ax  lies  to  the  left  of  the  rarefaction  wave,  then 

p  =  p^,  u  =  u^,  and  p  =  p^  o 
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(2)  If  i   Ax  lies  inside  the  right  rarefaction  wave,  we 

dx 
equate  the  slope  of  the  characteristic  -g^  =  ^  +  ^  to  the  slope  of 

the  line  through  the  origin  and  (I^^Ax,  At/2 ),  obtaining 

2i   Ax 
u+c=-^.  (57) 

With  the  constancy  of  the  Riemann  invariant 

2c(7-l)-l-  u  =  2c^(7-l)-^-u^  (58) 

the  isentropic  law,  cind  the  definition  of  c,  we  can  obtain  p,  u, 
and  p.   Using  the  isentropic  law  we  obtain 


P  =  PpPj,"^p^  =  Ap^  .  (59  ) 

Using  equation  (58 )  we  obtain,  by  solving  for  c 

c  =  c^+  ^^  (u  -  u^)  .  (60  ) 

Substitution  of  (60))  into  (.57)  and  solving  for  u  we  obtain 

o    2C  Ax        -, 

By  substitution  of  (g-i  0  into  (  60)   c  is  obtained;  by  substitution 
of  (  59)  into  the  definition  of  c  and  solving  for  p  we  obtain 

,c2lA-l 
P  =  (75^      *  (62) 

(3)  If  ^  Ax  lies  to  the  right  of  the  right  rarefaction  wave 
we  obtain  p  =  p  ,  u  =  u  ,  and  p  =  p  » 

Equations  {  i^c)   -    (  62)  are  the  key  to  the  programming  of 
Glimm's  method.   For  a  summary  see  the  flow  chart.  Fig. 16  • 
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Table  II 


X 

p 

u 

P 

e 

* 

1+  . 

0.1 

1.000 

0.000 

1.000 

2.500 

2.958 

0.2 

1.000 

0.000 

1.000 

2.500 

2.958 

0.3 

0.869 

0.164 

0.822 

2.363 

2.958 

0.^ 

0.'426 

0.927 

0.303 

1.778 

2.958 

0.5 

0.426 

0.927 

0.303 

1.778 

2.958 

0.6 

0.426 

0.927 

0.303 

1.778 

2.958 

0.7 

0.426 

0.927 

0.303 

1.778 

2.958 

0.8 

0.266 

0.927 

0.303 

2.853 

3.624 

0.9 

0.125 

0.000 

0.100 

2.000 

2.646 

*  r  is  the  Riemann  Invariant   ^31  ^  I  '  where  c  is  the  local 
sound  speed. 


32 


Table  III 


t/At 


Qr 


Q, 


m 


1 

Oo^^? 

O.OlB 

2.213 

2 

0.550 

0.019 

2.217 

3 

0.55'i 

0.032 

2.218 

i* 

0.550 

0.039 

2.219 

5 

0.552 

0.04? 

2.223 

6 

0.550 

0.059 

2.222 

7 

0.5^19 

0.070 

2.221 

8 

0.550 

0.079 

2.224 

9 

0.5^5 

0.090 

2.232 

10 

0.5^^6 

0.097 

2.247 

11 

0.5^18 

0.110 

2.258 

12 

0.5^5 

0.119 

2.267 

13 

0.5^9 

0.122 

2.266 

Ik 

0.552 

0.136 

2.266 

15 

0.5'<9 

0.1^13 

2.269 

16 

0.553 

0.1il9 

2.275 

17 

0.550 

0.158 

2.266 

18 

0.5'46 

0.164 

2.267 

19 

0.550 

0.178 

2.267 

20 

0.5^3 

0.190 

2.272 
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Table  IV 


SCHEMES 

WITHOUT  ACM 
0.226 

WITH  ACM 

Godunov 

0.2il7 

Lax-Wendroff 

0.226 

- 

MacCormack 

0.22i< 

- 

Rusanov 

0.22^4 

0.2^40 

Upwind 

0.225 

- 

Glimm 

0.36'4 

- 

Antidiffusion 

0.242 

- 

Hybrid 

0.258 

0.269 

Hyman 

0.276 

- 

Times  include  computation  of  exact  solution,  calls  to  printing  and 
plotting  routines,  which  were  the  same  for  all  cases. 
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Figure   6   continued 
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